#include "AnaDef.h"
#include "style.h"

const Double_t EPSILON = 1e-10;

void stylec0()
{
  style::fgkTitleSize = 0.07;
  style::fgkTextSize = 0.07;//for legend
  style::fgkYTitleOffset = 1;
}

void stylec1()
{
  style::fgkTitleSize = 0.1;//title
  style::fgkTextSize = 0.1;//label, legend (if any)
  style::fgkYTitleOffset = 0.7;
  style::fgkXTitleOffset = 1.35;
}

//3000
Double_t IDdEdxM3Q0(const Double_t *xx, const Double_t *par)
{
  const Double_t mean=xx[0];
  return                    TMath::Abs(par[0])/TMath::Sqrt(mean) + par[1] + par[2]/mean;
}

//2000
Double_t IDdEdxM2Q0(const Double_t *xx, const Double_t *par)
{
 const Double_t mean=xx[0];
  return                  ( TMath::Abs(par[0]) * TMath::Power(mean,par[1]) + par[2] )/mean;
}

Double_t RawAleph(const Double_t *xx, const Double_t *par)
{
  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = TMath::Abs(par[0]);
  const Double_t p1 = TMath::Abs(par[1]);
  const Double_t p2 = par[2];//TMath::Abs(par[2]);//still first have to be positive
  const Double_t p3 = TMath::Abs(par[3]);
  const Double_t p4 = TMath::Abs(par[4]);//has to be positive so that in bg-> inf there is no bg dependence anymore

  //return p0*((p1-TMath::Log(p2+1/TMath::Power(bg,p4)))/TMath::Power(beta,p3) - 1 );
  //return p0*( (p1- TMath::Log(p2+EPSILON)- TMath::Log(1+1/(TMath::Power(bg,p4)*p2+EPSILON))  )/( TMath::Power(beta,p3)+EPSILON ) - 1 );
  //log(p2) -> p2 because p2 sometimes ~EPSILON
  return p0*( (p1- p2 - TMath::Log(1+1/(TMath::Power(bg,p4)*TMath::Exp(p2)+EPSILON))  )/( TMath::Power(beta,p3)+EPSILON ) - 1 );
  //p2 can be very small ~ 1e-9, do not separate p2 out! -> absorb p2 into p1, remember to change initial value of p1
  //return p0*( (p1- TMath::Log(1+1/(TMath::Power(bg,p4)*p2+EPSILON))  )/( TMath::Power(beta,p3)+EPSILON ) - 1 );
  //change p0 as plateau
  //return p0/(p1-1)*( (p1- TMath::Log(1+1/(TMath::Power(bg,p4)*p2+EPSILON))  )/( TMath::Power(beta,p3)+EPSILON ) - 1 );
  //redefine p1 so that p1>1 --> p1>0
  //return p0/p1*( (p1 + 1 - TMath::Log(1+1/(TMath::Power(bg,p4)*p2+EPSILON))  )/( TMath::Power(beta,p3)+EPSILON ) - 1 );
}

Double_t myBetheBloch(const Double_t *xx, const Double_t *par)
{
  //RawALEPH has numerical difficulty in ln(p2+1/pow(bg, p4))
  //use flat dedx for electron. deviation from this at high bg is not inportant because anyway there the statistics is very limited
  //npar = 5;

  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = par[0];//overall prefactor
  const Double_t p1 = par[1];//beta exponent ~ 2
  const Double_t p2 = par[2];
  const Double_t p3 = par[3];
  const Double_t p4 = par[4];//plateau ~ 80

  //EPSILON added to protect 0 denominator!
  const Double_t bethebloch = p0 *( 1 + p2*TMath::Log(bg) + p3*TMath::Log(beta) )/( TMath::Power(beta, p1)+EPSILON) ;

  //tp = 100, electron momentum = 51 MeV, pion momentum = 14 GeV
  const Double_t tp = 10;

  Double_t qq = 0;
  if(bg<tp){
    qq = bethebloch;
  }
  else{
    qq = TMath::Min(p4, bethebloch);
  }

  return qq;
}

TF1 * GetmyBetheBloch(const TString tag, const Double_t pars[], const Int_t col, const Int_t sty, const Double_t xmin, const Double_t xmax)
{
  TF1 *f1=new TF1(Form("myBetheBloch%s", tag.Data()), myBetheBloch, xmin, xmax, 5);
  f1->SetLineColor(col);
  f1->SetLineStyle(sty);
  f1->SetParameters(pars);
  return f1;
}

TF1 * GetRawAleph(const TString tag, const Double_t pars[], const Int_t col, const Int_t sty, const Double_t xmin, const Double_t xmax)
{
  TF1 *f1=new TF1(Form("RawAleph%s", tag.Data()), RawAleph, xmin, xmax, 5);
  f1->SetLineColor(col);
  f1->SetLineStyle(sty);
  f1->SetParameters(pars);
  return f1;
}

TF1 * GetRes3000(const TString tag, const Double_t pars[], const Int_t col, const Int_t sty, const Double_t xmin, const Double_t xmax)
{
  TF1 *f1=new TF1(Form("Res3000%s", tag.Data()), IDdEdxM3Q0, xmin, xmax, 3);
  f1->SetLineColor(col);
  f1->SetLineStyle(sty);
  f1->SetParameters(pars);
  return f1;
}

TF1 * GetRes2000(const TString tag, const Double_t pars[], const Int_t col, const Int_t sty, const Double_t xmin, const Double_t xmax)
{
  TF1 *f1=new TF1(Form("Res2000%s", tag.Data()), IDdEdxM2Q0, xmin, xmax, 3);
  f1->SetLineColor(col);
  f1->SetLineStyle(sty);
  f1->SetParameters(pars);
  return f1;
}

void drawFunc(TPad * pad0, TPad * pad1, TF1 *ff, TF1 *fref, const TString dopt, const TString ytit, const Double_t ymin, const Double_t ymax, const Double_t rmin, const Double_t rmax)
{
  stylec0();
  style::PadSetup(pad0);
  pad0->SetLeftMargin(0.14);
  pad0->SetRightMargin(0.03);
  pad0->SetTopMargin(0.1);
  pad0->SetBottomMargin(0.01);
  pad0->cd();

  //histogram style works normal
  TH1D * h0 = style::Func2Hist(ff, kTRUE);
  style::ResetStyle(h0);
  h0->SetMinimum(ymin);
  h0->SetMaximum(ymax);
  h0->SetLineWidth(2);
  h0->Draw(dopt);

 //----------------------
  stylec1();
  style::PadSetup(pad1);
  pad1->SetTopMargin(0.01);
  pad1->SetLeftMargin(0.14);
  pad1->SetRightMargin(0.03);
  pad1->SetBottomMargin(0.3);
  pad1->cd();

  //GetHistogram very very bad
  TH1D * hf = style::Func2Hist(ff, kTRUE); //(TH1D*)ff->GetHistogram()->Clone(Form("%sff",ff->GetName()));
  TH1D * href = style::Func2Hist(fref, kTRUE);//fref->GetHistogram();

  hf->Divide(href);
  style::ResetStyle(hf);

  hf->SetMinimum(rmin);
  hf->SetMaximum(rmax);
  hf->GetXaxis()->SetTitle(ff->GetXaxis()->GetTitle());
  hf->GetYaxis()->SetTitle(ytit);
  hf->GetYaxis()->SetNdivisions(505);
  hf->SetLineWidth(2);
  hf->Draw(dopt);

}

void drawdEdxTPCCutMode01(const Int_t kMC=0)
{
  style::SetGlobalStyle(0);

  TCanvas *cc=new TCanvas("cc","",600,500);

  TPad *c0=new TPad("c0","",0,0.4,1,1);
  c0->SetLogx(1);
  c0->Draw();

  TPad *c1= new TPad("c1","", 0,0,1,0.4);
  c1->SetLogx(1);  
  c1->Draw();

  //-----------------------------------
  Double_t fPars[10];

  //all checked no bad fit

  //DATA
  //xlu@lxsub01:/hera/alice/xlu/workspace/hera2dbetabloch/FracFunc_10-d-e/more_datajet_10032014/f09finish$ for ii in $(echo TPC*TOF9/*ETA0*JPT5-*/XB5YB2*mBB*anchor3*RES3*/coherentFit/logmn2d.txt ); do echo $ii ; cat $ii | grep -vi bin | tail -n 21 | head -n 11;done
  /*
TPC1TOF9/DATA_histRecTHnIDFF_LHC10de_TPC1TOF9_kMC0_ETA0_JPT5-10_coherentg09/XB5YB2_Xmin-9_Xmax40_mBB_anchor3_NTYPE4_SignalShift_RES3000_POLINT2_Tol2.0_RegOpt0__FixHighPElectronMode2/coherentFit/logmn2d.txt
fPars[0] = 3.870488e+01;
fPars[1] = 2.077405e+00;
fPars[2] = 1.590383e-01;
fPars[3] = -2.296176e-02;
fPars[4] = 7.386469e+01;
fPars[5] = 0.000000e+00;
fPars[6] = 0.000000e+00;
fPars[7] = -3.592465e-03;
fPars[8] = 5.792973e-01;
fPars[9] = -1.365681e+01;
fPars[10] = 5.150183e+01;
TPC2TOF9/DATA_histRecTHnIDFF_LHC10de_TPC2TOF9_kMC0_ETA0_JPT5-10_coherentf09/XB5YB2_Xmin-9_Xmax40_mBB_anchor3_NTYPE4_SignalShift_RES3000_POLINT2_Tol2.0_RegOpt0__FixHighPElectronMode2/coherentFit/logmn2d.txt
fPars[0] = 3.857770e+01;
fPars[1] = 2.148615e+00;
fPars[2] = 1.603453e-01;
fPars[3] = 3.960138e-02;
fPars[4] = 7.390617e+01;
fPars[5] = 0.000000e+00;
fPars[6] = 0.000000e+00;
fPars[7] = -3.565727e-03;
fPars[8] = 5.709401e-01;
fPars[9] = -1.363173e+01;
fPars[10] = 5.119855e+01;

  */
  //MC
  //
  /*
TPC1TOF9/DATA_histRecTHnIDFF_LHCMC10de_TPC1TOF9_kMC1_ETA0_JPT5-10_coherentg09/XB5YB2_Xmin-9_Xmax40_mBB_anchor3_NTYPE4_SignalShift_RES3000_POLINT2_Tol2.0_RegOpt0__FixHighPElectronMode2/coherentFit/logmn2d.txt
fPars[0] = 3.418394e+01;
fPars[1] = 2.427050e+00;
fPars[2] = 1.929519e-01;
fPars[3] = 3.182934e-01;
fPars[4] = 7.339055e+01;
fPars[5] = 0.000000e+00;
fPars[6] = 0.000000e+00;
fPars[7] = -1.040427e-02;
fPars[8] = 3.450450e-01;
fPars[9] = -1.291983e+01;
fPars[10] = 4.918995e+01;
TPC2TOF9/DATA_histRecTHnIDFF_LHCMC10de_TPC2TOF9_kMC1_ETA0_JPT5-10_coherentf09/XB5YB2_Xmin-9_Xmax40_mBB_anchor3_NTYPE4_SignalShift_RES3000_POLINT2_Tol2.0_RegOpt0__FixHighPElectronMode2/coherentFit/logmn2d.txt
fPars[0] = 3.411760e+01;
fPars[1] = 2.452716e+00;
fPars[2] = 1.939329e-01;
fPars[3] = 3.347598e-01;
fPars[4] = 7.344232e+01;
fPars[5] = 0.000000e+00;
fPars[6] = 0.000000e+00;
fPars[7] = -1.066311e-02;
fPars[8] = 2.745010e-01;
fPars[9] = -1.223212e+01;
fPars[10] = 4.590847e+01;

*/

  //-----------------
  if(kMC==0){
    fPars[0] = 3.870488e+01;                                                                                                                                                                                                                                     
    fPars[1] = 2.077405e+00;                                                                                                                                                                                                                                     
    fPars[2] = 1.590383e-01;                                                                                                                                                                                                                                     
    fPars[3] = -2.296176e-02;                                                                                                                                                                                                                                    
    fPars[4] = 7.386469e+01;                                                                                                                                                                                                                                     
    fPars[5] = 0.000000e+00;                                                                                                                                                                                                                                     
    fPars[6] = 0.000000e+00;                                                                                                                                                                                                                                     
    fPars[7] = -3.592465e-03;                                                                                                                                                                                                                                    
    fPars[8] = 5.792973e-01;                                                                                                                                                                                                                                     
    fPars[9] = -1.365681e+01;                                                                                                                                                                                                                                    
    fPars[10] = 5.150183e+01;     
  }
  else{
    fPars[0] = 3.418394e+01;                                                                                                                                                                                                                                     
    fPars[1] = 2.427050e+00;                                                                                                                                                                                                                                     
    fPars[2] = 1.929519e-01;                                                                                                                                                                                                                                     
    fPars[3] = 3.182934e-01;                                                                                                                                                                                                                                     
    fPars[4] = 7.339055e+01;                                                                                                                                                                                                                                     
    fPars[5] = 0.000000e+00;                                                                                                                                                                                                                                     
    fPars[6] = 0.000000e+00;                                                                                                                                                                                                                                     
    fPars[7] = -1.040427e-02;                                                                                                                                                                                                                                    
    fPars[8] = 3.450450e-01;                                                                                                                                                                                                                                     
    fPars[9] = -1.291983e+01;                                                                                                                                                                                                                                    
    fPars[10] = 4.918995e+01;      
  }

  TF1 * fmeanTPC1=GetmyBetheBloch("meanTPC1", fPars, kRed, kSolid, 0.4, 1e4);
  fmeanTPC1->GetXaxis()->SetTitle("#beta#gamma");
  fmeanTPC1->GetYaxis()->SetTitle("#mu (arb. units)");

  TF1 * fsigmaTPC1=GetRes3000("sigmaTPC1", &(fPars[8]), kRed, kSolid, 40, 90);
  fsigmaTPC1->GetXaxis()->SetTitle("#mu (arb. units)");
  fsigmaTPC1->GetYaxis()->SetTitle("#sigma'");

  //-----------------
  if(kMC==0){
    fPars[0] = 3.857770e+01;                                                                                                                                                                                                                                     
    fPars[1] = 2.148615e+00;                                                                                                                                                                                                                                     
    fPars[2] = 1.603453e-01;                                                                                                                                                                                                                                     
    fPars[3] = 3.960138e-02;                                                                                                                                                                                                                                     
    fPars[4] = 7.390617e+01;                                                                                                                                                                                                                                     
    fPars[5] = 0.000000e+00;                                                                                                                                                                                                                                     
    fPars[6] = 0.000000e+00;                                                                                                                                                                                                                                     
    fPars[7] = -3.565727e-03;                                                                                                                                                                                                                                    
    fPars[8] = 5.709401e-01;                                                                                                                                                                                                                                     
    fPars[9] = -1.363173e+01;                                                                                                                                                                                                                                    
    fPars[10] = 5.119855e+01;  
  }
  else{
    fPars[0] = 3.411760e+01;                                                                                                                                                                                                                                     
    fPars[1] = 2.452716e+00;                                                                                                                                                                                                                                     
    fPars[2] = 1.939329e-01;                                                                                                                                                                                                                                     
    fPars[3] = 3.347598e-01;                                                                                                                                                                                                                                     
    fPars[4] = 7.344232e+01;                                                                                                                                                                                                                                     
    fPars[5] = 0.000000e+00;                                                                                                                                                                                                                                     
    fPars[6] = 0.000000e+00;                                                                                                                                                                                                                                     
    fPars[7] = -1.066311e-02;                                                                                                                                                                                                                                    
    fPars[8] = 2.745010e-01;                                                                                                                                                                                                                                     
    fPars[9] = -1.223212e+01;                                                                                                                                                                                                                                    
    fPars[10] = 4.590847e+01;  
  }

  TF1 * fmeanTPC2=GetmyBetheBloch("meanTPC2", fPars, kBlue, kDotted, 0.4, 1e4);
  fmeanTPC2->GetXaxis()->SetTitle("#beta#gamma");
  fmeanTPC2->GetYaxis()->SetTitle("#mu (arb. units)");

  TF1 * fsigmaTPC2=GetRes3000("sigmaTPC2", &(fPars[8]), kBlue, kDotted, 40, 90);
  fsigmaTPC2->GetXaxis()->SetTitle("#mu (arb. units)");
  fsigmaTPC2->GetYaxis()->SetTitle("#sigma'");
  
  //========================================================

  TLegend *lg=new TLegend(0.55, 0.65, 0.98, 0.88);
  stylec0();
  style::ResetStyle(lg);

  //drawFunc(c0,c1, fmeanTPC0, fmeanTPC0, "l", "ratio to ref.",              45, 135, 0.995, 1.005);
  drawFunc(c0,c1, fmeanTPC1, fmeanTPC1, "l", "ratio to ref.",      45, 135, 0.995, 1.005);
  drawFunc(c0,c1, fmeanTPC2, fmeanTPC1, "l same", "",      45, 135, 0.995, 1.005);

  c0->cd();

  lg->SetHeader(Form("%s jet pt 5-10 GeV/c tracks, 0.4<|#eta|<0.6", kMC? "MC 10d":" data 10d+e"));
  //lg->AddEntry(fmeanTPC0,"hybrid cut (ref.)","l");
  lg->AddEntry(fmeanTPC1,"TPCsignalN>60 (ref.)", "l");
  lg->AddEntry(fmeanTPC2,"TPC Geo.", "l");
  lg->Draw();

  cc->Print(Form("TPCcutModemeankMC%d.eps", kMC));
  //========================================================
  lg=new TLegend(0.55, 0.65, 0.98, 0.88);
  stylec0();
  style::ResetStyle(lg);

  //drawFunc(c0,c1, fsigmaTPC0, fsigmaTPC0,   "l", "ratio to ref.",   0.039, 0.081, 0.88, 1.12);
  drawFunc(c0,c1, fsigmaTPC1, fsigmaTPC1,   "l", "ratio to ref.",           0.039, 0.081, 0.88, 1.12);
  drawFunc(c0,c1, fsigmaTPC2, fsigmaTPC1,   "l same", "",           0.039, 0.081, 0.88, 1.12);
 
  //c0->SetTopMargin(0.08);
  c0->SetLogx(0);
  c1->SetLogx(0);
  c0->cd();

  lg->SetHeader(Form("%s jet pt 5-10 GeV/c tracks, 0.4<|#eta|<0.6", kMC? "MC 10d":" data 10d+e"));
  //lg->AddEntry(fsigmaTPC0,"hybrid cut (ref.)","l");
  lg->AddEntry(fsigmaTPC1,"TPCsignalN>60 (ref.)", "l");
  lg->AddEntry(fsigmaTPC2,"TPC Geo.", "l");
  lg->Draw();

  cc->Print(Form("TPCcutModesigmakMC%d.eps", kMC));
}

int main()
{
  drawdEdxTPCCutMode01();
  return 0;
}
